#include "GK_PC_solver.h"

/* ************************************************************************************************  *
                                         GKPC Solver 
*  ************************************************************************************************* */

template<class T>
long countNZ(T *ptr, size_t size)
{
	long nz = 0;
	for(index i=0;i<size;i++)
			nz += ptr[i]!=0 ? 1 : 0;
	
	
	return nz;
}





void calculateP(matrix_type *A,real *b,dreal *x, dreal *y, dreal *P,dreal delta, index m, index n)
{
	//calculates P[q] = min_i b(i)/A(i,q), q=1,2,...,m  
	//                -> matrix A in a column-wise representation
	
    dreal min;
	dreal val;
	
	// initialize arrays x,y
	for(index i=0;i<n;i++) x[i] = 0.0;
	for(index i=0;i<m;i++) 
		if (b[i]>0) y[i] = delta/b[i];
		else y[i] = 0;
	
	for(index col=0;col<n;col++)
	{
	   min = INFTY;
	   for(index row=0;row<m;row++)
	   {
	      val=abs(A[row+col*m])>0 ? b[row]/abs(A[row+col*m]):INFTY;
	      if(val<min) min = val;
		}
	P[col] = min;
		
    }

};

void findRho(matrix_type *A,real *c,dreal *y,dreal &rho, index &q, integer m, integer n)
{
	//calculates R[j] = sum_i A[i][j] y[i] /c[j], j=1,2,...,m i.e. length_y(j)
	rho = INFTY;
	q = 0;
    
	dreal val;
	for(index col=0;col<n;col++)
	{
		val = 0.0;
	   for(index row=0;row<m;row++)
	   {
	      val+=abs(A[row+col*m])*y[row];
	   }
	   val /= c[col];
	   
	   if(val<rho && val>0)
	   {
			rho = val;
			q   = col;
			
	   }
	   
    }
	
};



void updatePrimalDual(matrix_type *A, real *b,real *c,dreal *x,dreal *y, dreal *P,dreal &primal, dreal &dual,dreal eps, index q,dreal rho, integer m)
{
   x[q] += P[q];
   primal += c[q]*P[q];
   dual += eps*c[q]*P[q]*rho;
   
   for(index i=0;i<m;i++)
   
   {
	if(y[i]!=0)
		y[i]*=(1.0+eps*P[q]*abs(A[i+q*m])/b[i]); 
	else y[i]=0;
   }
}

void scalePrimalDual(dreal *x, dreal *y, dreal primalScale, dreal dualScale,integer m, integer n)
{
     for(index i=0;i<n;i++) x[i] /= primalScale;
     for(index i=0;i<m;i++) y[i] /= dualScale;
}





LP_status GKPCSol(matrix_type *A, real *b, real *c, dreal *x, dreal *y, integer m, integer n,dreal omega)
{
	if(VERBOSE)
	{
	std::cout << std::endl << std::endl;
	std::cout << "---------------------------------------------------------------" << std::endl;
	std::cout << "         Garg-Koenemann packing-covering LP solver               " << std::endl;
	std::cout << "---------------------------------------------------------------" << std::endl;
	}

	// count non-zero elements
	long nz = countNZ(A,m*n);
	
	Timer compute_time;
	
	compute_time.start();
	
	// parameters defined by algorithm
	//dreal eps = 1-pow(1.0+omega, (double) -1.0/1.6);
	dreal eps = 3.0/2.0 - sqrt(1.0+8.0/(omega+1.0))/2.0;		
	dreal delta = (1.0+eps)*pow(((1+eps)*m),(-1.0/eps));
	
	
	// set initial primal dual value
	dreal primal = 0.0;
        dreal dual   = 0.0;
      
	for(index i=0;i<m;i++) if(b[i]>0) dual+=delta/b[i]; // exchange for dual = m*delta
	
	// scale factors
	dreal primalScale = 1.0-(log(delta)/log(1+eps));
	dreal dualScale = 1.0; // =(*rho) later

	// 
	dreal *P = new dreal[n];
	dreal rho = 0;
	index q;
	
	if(VERBOSE)
	{
	std::cout << "---------------------- parameters -----------------------------" << std::endl;
	std::cout << "omega: " << omega << "\t eps: " << eps << "\t delta: " << delta << std::endl; 
	}
	
	// PREPROCESSING
	if(VERBOSE) std::cout << "preprocessing step ..." << std::endl;
	calculateP(A,b,x,y,P,delta,m,n);
	
	if(VERBOSE) std::cout << "done." << std::endl; 
	integer iter = 0;
	
	
	// iterate to find primal-dual pair (x,y)
	if(VERBOSE) std::cout << "finding approximate solution ..." << std::endl;
	
	while(dual<1.0)
	{
		iter ++;
		// find rho
		findRho(A,c,y,rho,q,m,n);
				
		// update primal-dual value
		updatePrimalDual(A,b,c,x,y,P,primal,dual,eps,q,rho,m);
			
	
	} // end_while
	
	
	if (VERBOSE) std::cout << "done." << std::endl;
	findRho(A,c,y,rho,q,m,n);
	if (VERBOSE) std::cout << "scaling dual solution with rho=" << rho << std::endl;
	dualScale = rho;

	// scaling final solution
	scalePrimalDual(x,y,primalScale,dualScale,m,n);
	primal /= primalScale;
	dual   /= dualScale;
		
	if(VERBOSE) std::cout << "done." << std::endl;
	
	compute_time.stop();

    if(VERBOSE)
	{
    std::cout << "\nDATA SET: rows: " << m << ", cols: " << m << ", density: " << nz/(double)(n*m)*100<< "%" << std::endl;
	std::cout << "---------------------------------------------------------------" << std::endl;
			
	std::cout << "Number of iterations:                           | " << iter << std::endl;
	std::cout << "----------------------------------------------------------------" << std::endl;
	std::cout << "*** Compute time:                               | " << compute_time.secs() << std::endl;
	std::cout << "----------------------------------------------------------------" << std::endl;

	printf(" primal value: %f       \n\r dual value: %f     \n\r approx. ratio: %f \n",primal,dual,dual/primal);
	std::cout << "----------------------------------------------------------------" << std::endl;
	}

	delete [] P;
	
	if(dual/primal<=1+omega) return LP_FEAS;
	else return LP_INF;
	
}
